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I. INTRODUCTION 

When dealing with measurements on complex systems it 
is typically difficult to find the optimal set of variables to de- 
scribe their evolution, which provides the best opportunities 
for understanding and predicting the system's behavior 

Recently a framework for analyzing measurements on 
stochastic systems was introduced Jll |21 ■ This framework is 
based on the assumption that the measured properties evolve 
according to a deterministic drift and stochastic fluctuations 
due to the interactions with the internal degrees of freedom 
or the environment. These fluctuations are globally ruled by 
some specific stochastic forcing inherent to the signal evolu- 
tion itself, and therefore cannot be separated from the mea- 
sured variables. Whether they are due to quantum uncertain- 
ties or a limited knowledge of the state of the system, they are 
fundamental for a complete description of a complex system. 
Several works in this scope have already shown the advan- 
tage of this approach, ranging from the description of turbu- 
lent flowsjfll] and climate indices [jstl to the evolution of stock 
markets jjand oil prices [jSI, just to mention some. In the 
course of time several improvements to the method were pro- 
posed concerning its robustness with respect to finite sampling 
effects and measurement noise ll6l- fTTIl . 

However, a natural question arises when dealing simultane- 
ously with several variables: is it possible to find non-trivial 
functions of measured properties for which the stochastic fluc- 
tuation can be neglected? This would mean that it should be 
possible to decrease the number of stochastic variables needed 
to describe the system. 

In this paper we will follow this question up in detail and 
show, that by accessing the eigenvalues of the diffusion ma- 
trices comprehending the measured properties it is possible to 
derive a path in phase space through which the deterministic 
contribution is enhanced. As a direct application, our proce- 
dure allows for the determination of the number of indepen- 
dent sources of stochastic forcing in a system described by an 



arbitrarily large number of properties. 



A pictorial example is as follows. Consider a professional 
archer trying to aim at a target in a succession of shots. If 
the archer holds the bow without any support, one expects 
that the set of trials is distributed around the center according 
to a radially symmetric Gaussian. Hence deviations have the 
same amplitudes in all directions. However, if the archer lies 
on the floor, the vertical direction will be more confined than 
the horizontal. In this case we expect a smaller variance of 
deviations in the vertical direction than in the horizontal one. 
If the archer lies on an inclined plane the most confined and 
less confined directions will have a certain slope related to 
the plane inclination. Regarding the variance in space as a 
particular representation of the diffusion matrix, we can now 
think in more complex systems where at each point in phase 
space fluctuations can be defined through their local diffusion 
matrices. The study of its eigenvalues and eigenvectors is the 
scope of the present paper. 



We start in Sec. |ll] by defining our system of variables 
mathematically and by describing the standard approach to 
estimate its dynamics from measured data as introduced in 
Ref. 11101, where a particular emphasis will be given to the 
diffusive terms. Then the general eigenvalue problem is intro- 
duced in the role of the local eigenvalues of diffusion matrices 
are discussed, which exhibits the framework for the present 
work. We show that at each point of phase space the eigen- 
vectors of the diffusion matrix are tangent to new coordinate 
lines, one of them corresponding to the lowest eigenvalue, 
thus indicating the direction towards which fluctuations are 
- at least partially - suppressed. In Sec. |III]we demonstrate 
our approach on two examples from Hopf-bifurcation systems 
with stochastic forcing. Section|IV]cIoses the paper with dis- 
cussion and conclusions. 
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II. DEFINING STOCHASTIC EIGENDIRECTIONS 

We consider an A^-dimensional Langevin process X = 
{Xi{t), . . . , X]\[{t)) whose probability density functions 
(PDFs) /(X, t) evolve according to the Fokker- Planck equa- 
tion (FPE) 



9/(X,t) 
dt 
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The functions and D^j' are called the Kramers-Moyal or 
the drift and diffusion coefficients and are defined as 

D«(X) = l.„, i-""''^'^" , ,2) 

^ ' At^o At kl 

where M^'^) are the first and second conditional moments 
{k — 1, 2). Here we assume that the underlying process is 
stationary and therefore both drift and diffusion coefficients 
do not explicitly depend on time t. Conditional moments can 
be directly derived from the measured data as ^ 



dxidxj 

-,(2) 



i?(^'(X)/(X,i) 



(1) 



= {Y,{t + M)-Y,{t)\Y{t)^lC) 

y,(t))(y,(t + At)-r,(t))|Y(f) = x) 



Mf^(X,At) 
(X,At) 
((F,(t + At)- 

(3) 

Here Y{t) = {Yi{t),..., YnII:)) exhibits the A^-dimensional 
vector of measured variables at the t and (•|Y(t) = X) sym- 
bolizes a conditional averaging over the entire measurement 
period, where only measurements with Y(t) = X are taken 
into account ll27ll . D'^' is the drift vector and D*^^' the diffu- 
sion matrix. 

Associated with the FPE ([T]) is a system of N coupled Ito- 
Langevin equations, which can be written as lfT2l[T3ll 



G(x)r(0 



(4) 



Here T{t) is a set of N normally distributed random variables 
fulfilling 



(r,(t)) = 0, {T,{t)V,{t')) = 2%<5(t - t') 



(5) 



that drives the stochastic evolution of X. The vectors h and 
the matrices G ~ {gij} for alH, = 1, . . . , are connected 
to the local drift and diffusion function through 



DY'{X) = h,{X) and 

N 

i^g'(X) = ^g,fe(X)g,,(X) 

k=l 



(6) 
(7) 



While the FPE describes the evolution of the joint distribu- 
tion of the N variables statistically, the system of Langevin 
equations in (|4]i models individual stochastic trajectories of 
the system. In Eq. (|4| the term h(X) contains the determinis- 
tic part of the macroscopic dynamics, while the locally linear 



functions G(X) account for the amplitudes of the stochastic 
forces mirroring the different sources of fluctuations due to all 
sorts of microscopic interactions within the system. Here we 
should point out that we consider stationary processes, other- 
wise ensemble averages have to be taken L2i|. 

The conditional moments in Eq. (O are computed directly 
from data time series and are typically linear functions of At 
for sufficiently small At 12|,|7l]. Then, through the limit in 
Eq. both Lif ^ and D^f are determined. 

The N X N matrix G cannot be uniquely determined from 
the symmetric diffusion matrix D'^' for N > 2, because the 
number of unknown elements in G exceeds the number of 
known elements in D^^' leading to A^^ - ^N{N + 1) = 
^N{N — 1) free parameters. However, a simple method to 
obtain G from D^^^ is the following. Due to its symmetry 
and positive semi-definiteness the diffusion matrix D^^^ has 
only real, non-negative eigenvalues Xi{i = 1, . . . , N). There- 
fore an orthogonal transformation U can be found that diag- 
onalizes D^^), i.e. U'^Dt^^U = diag(Ai, . . . , An). Taking 
the positive root of the eigenvalues and transforming it back 
we arrive at gij = (vD^)ij where the symbolic notation 

)U-^ is used for simplicity. 



VdW = Udiag(VA7, ...,V^N) 
General forms of G can be constructed by multiplication of 
\/D(2) with arbitrary orthogonal matrices. For our consider- 
ations the simple version is sufficient lfl2ll . 




FIG. 1: Stream plot of the system in Eq. (|8} with no stochasticity 
(D'^' = 0), i.e. kl = k'z = 0, where the aiTows formally represent 
the drift vector of synthetic data extracted from direct simulation of 
Eq. {8](. The solid line represents a partial trajectory. Trajectories 
diverge from the unstable fixed point r — and converge to the limit 
cycle r = 1. 
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FIG. 3: Time series of (a) x{t) = r{t) cos e{t) and (b) y{t) = 
r{t) sin 9{t) where r and ^ are talcen from the integration of Eq. (|8) 
for the same conditions as in Fig.|2] 



more, if the corresponding eigenvalues are very small com- 
pared to all the other ones at the respective mesh points, the 
corresponding stochastic forces can be neglected and the sys- 
tem has only iV — 1 independent stochastic forces. In this case 
the problem can be reduced in one stochastic variable by an 
appropriate transformation of variables, since the eigenvectors 
in one coordinate system are the same as in another one (see 
Append, [a]). 

These are the central ideas of our study which we next apply 
to an analytical example, namely the Hopf -bifurcation. 



FIG. 2: Principal axis of the diffusion matrix. Field of (a) largest and 
(b) smallest eigenvectors for the diffusion matrix corresponding to 
Eq. ^ together with (c) the eigenvalues A as a function of r and (in- 
set) a close-up around A ~ (see text). The largest eigenvalue is as- 
sociated with the radial direction while the smallest one corresponds 
to the angular one. Here fci = 0.5, k2 = 0.05, a = 0.7475 ^ (r^), 
to = 0, ro = 1.01 and 6*0 ~ 0. The analysis was performed from 
10^ data points extracted by integration of Eq. ([8) with integration 
step of At = 10"''. 



Since h and G are functions of the variables X and are 
numerically determined onarii x ... x njv mesh of points 
in phase space, one can always define at each mesh point the 
N eigenvalues and corresponding eigenvectors of the matrix 
G(X). This analysis provides information about the stochas- 
tic forcing acting on the system and was already applied to a 
two-dimensional sub-critical bifurcation ifisll and to the anal- 
ysis of human movements lfl6ll . 

In general the eigenvalues indicate the amplitude of the 
stochastic force and the corresponding eigenvector indicates 
the direction toward which such force acts. Even more inter- 
esting features, however, can be extracted from the eigenval- 
ues and eigenvectors. 

To each eigenvector of the diffusion matrix we can asso- 
ciate one independent source of stochastic forcing F;. In this 
scope, the eigenvectors can be regarded as defining principal 
axes for stochastic dynamics. For instance, the vector field 
aligned at each mesh point to the eigenvector associated to 
the smallest eigenvalue of matrix G defines the paths in phase 
space towards which the fluctuations are minimal. Further- 



Ill. THE HOPF-BIFURCATION SYSTEM WITH 
STOCHASTIC FORCING: AN ANALYTICAL EXAMPLE 

In what follows we consider the dynamical system 

(IT 

= r(l - r2) + fcirFi (8a) 

at 

^ = a-r^ + fcsFs (8b) 
at 

describing the evolution of the radial and azimuthal coor- 
dinates of a particle moving in two-dimensional space, with 
a, ki and k2 being constants, ri.2 are the stochastic forces, 
which are Gaussian distributed and ^-correlated. For fci = 
k2 — the system for a > has an unstable fixed point at 
r = and a stable limit cycle at r = 1, as sketched in Fig.[T] 
Equation (|8bl ) uses a — instead of the usual 1 + r'^ lUTIl in 
order to be able to avoid the systematic increase of 6 in time. 
For this reason we chose a ~ (r^). 

Introducing non-zero stochastic terms ki ,k2^0 for inde- 
pendent stochastic forces, Fi and one has diagonal dif- 
fusion matrices with two eigenvalues, (fcir)^ and fc| at each 
point (r, 6) corresponding to the radial and angular eigendi- 
rections respectively. 

Since the eigenvectors do not change when changing the 
coordinate system (see Append. |A]i, "mixing" the radial and 
angular coordinates by e.g. changing to Cartesian should yield 
the same principal axis. We, therefore, now repeat our analy- 
sis for transformed variables x = r cos 9 and y ^ r sin 9. 

Using Ito's formulation for the transformation from polar 
to Cartesian coordinates [18:], the Langevin system can be 





FIG. 4: Drift and diffusion functions for tlie system in Eq. ^ for tfie same conditions as in Fig.|2l obtained from analysis of tiie time series 
X = rcos9 and y = rsinO. For eacfi drift and diffusion function tlie bottom surface reflect the numerical results, while the corresponding 
fitted surface on top is shifted upwards for clarity: (a) {x, y), (b) {x, y), (c) {x, y), (d) D[^^ {x, y), (e) -Dj^' (x, y). The bar plots 
compare the coefficients that fit each function (dark bars) with the analytic, true ones (white bars) (see Tab.|l). 



defined through the drift and diffusion coefficients in trans- 
formed coordinates, hi{x, y),gij{x, y) (see Append. lAb: 



= {h,/r~g|^)Y, + i-^yh2iYl^2^ + Y2^u)(9a) 

~g,j = ^{YiSa + Y2d,2) - g2J{Yl^2^ + Y2Su) (9b) 

with i,j = 1,2 and (Yi,F2) = {x,y). Thus, from G in 
the polar coordinate system, Eq. (|9Fi yields G for which the 



diffusion matrix D'^^ = GG reads 



D 



(2) 



-klxy - 



- klxy 



-klxy 



\- kfxy 



(10) 



The eigenvalues and corresponding eigenvectors are Ai = 
j/^) with vi = [;xl \fxP^~\^ ,yl \fxP~^r]p) and 
j(a;^ +t/^) with V2 



klix' 

A2 



In contrast to the eigendirections, the eigenvalues depend 
on the Jacobian of our transformation, since the nonlinear 
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TABLE I: Coefficients of quadratic forms describing the surfaces in Fig.|4] 



transformation from polar to Cartesian coordinates changes 
the metric. According to Append. |A] the eigenvalues in po- 
lar coordinates are = / sf for i = 1,2 with = 
(f )2 + ^ 1 and si = (If )2 + (If )2 = r\ The 

eigenvalues of our system in Cartesian coordinates are shown 
as solid and dashed lines in Fig. |2] and a full derivation for the 
eigenvalues in different coordinate systems is given in Ap- 
pend. |A]for the general iV-dimensional case. 

Notice that, if we determine the eigenvalues in a Cartesian 
system, the metric is Euclidean, i.e. the eigenvalues measured 
in the same units in both x- and y-direction directly charac- 
terize the diffusion in principal directions. Nonlinear transfor- 
mations of coordinates such as the one from Cartesian to polar 
coordinates, however, generally change the metric. In such a 
system the direction of the maximal eigenvalue is not neces- 
sarily the direction with the highest diffusion. For example, in 
the Hopf-bifurcation in Eqs. (O, the maximal eigenvalue for 
r < A:2 /fci is in azimuthal direction, but the maximal diffusion 
is still in the radial direction. 

Looking to x and y as measured variables, they define our 
Cartesian system and therefore Ai is our maximal eigenvalue 
and A2 is associated with the eigendirection showing minimal 
stochastic fluctuation. 

We proceed to verify this result by numerically analyzing 
the time series x and y plotted in Fig. [3j following the pro- 
cedure described in the previous section. The procedure is 
as follows. First, we compute the time series according to 
Eq. ^ in coordinates r and 9. Then, we transfer this time se- 
ries to Cartesian coordinates x, y. In the new coordinates we 
calculate the drift vectors D'^^ and diffusion matrices D*^^'. 
Finally, having matrices D*-^-* at each mesh point, one eas- 
ily computes the eigenvalues and eigenvectors, as shown in 
Fig.|2] The analytical results given by Eq. ( fTOl i are reproduced 
nicely: the stochastic contribution has two eigendirections, 
one in the radial direction and another in the angular one; the 
larger eigenvalue is associated with the radial direction. The 
deviations observed for the minimum eigenvalue occur due to 
differences in the quality of statistics in different regions of 
phase space (value of r): in particular regions in phase space 
at intermediate values of r are less frequently realized in the 
data series than others, resulting in a reduced accuracy of the 
estimated eigenvalues in these regions. 

Figures |4^-c show the components of the diffusion matrix, 

~ (2) ~ (2) ~ (2) 

namely D^^ , and D22 as a function of both x and y. 
Notice that the diffusion matrix is symmetric (see Eq. (|3)) and 



therefore D. 



(2) 
21 



■ The surfaces obtained numerically 



are fitted to a full quadratic polynomial in x and y through a 
least square procedure. The coefficients of these polynomial 
fits are also shown in Fig|4^-c. The results of Fig.|4]i-c comply 
with the analytical expression in Eq. (fTOl i. A similar analysis 
is done for the drift vector, i.e. for the function d[^'^ and D2^\ 
shown in Figs.HJl-e. Table U lists the values obtained for the 
coefficients of the surface fits in Fig.|4] 

It is worth stressing that the aforementioned procedure can 
likewise be applied to measured multi-dimensional time series 
in cases where no additional information on the underlying 
dynamics is available. In this case the corresponding fields of 
eigenvectors indicate the directions with largest and smallest 
stochastic contributions. Using this information, a transform 
of variables to a coordinate system aligned with the direction 
of the smallest stochastic contribution can be applied, which 
in the present case would be the tangential direction. In the 
representation of these new coordinates, the stochastic contri- 
bution will then be reduced. 

In the previous case, the two eigenvectors for diffusion were 
introduced tangentially and perpendicularly to the limit cy- 
cle. If the principal axes for the stochastic contribution are not 
aligned with the trajectories of the deterministic part of the 
system, one might suspect an interplay between both stochas- 
tic contribution and the drift of the system. 

In order to demonstrate the wider applicability of our 
method, we next consider such a case, with the major axis 
aligned with y-direction. We therefore investigate the system 



dr 
dd 



= r{l-r^) + Kircos9Ti+K2riiineT2 (11a) 
= (a -r2)-A'isin6'ri+A'2 008 6^2 .(lib) 



In Cartesian coordinates the dynamics are described by 

dx 



— = hi+KirVi 
at 

^ = /i2 + A'2rr2. 
dt 



with 



+ m - 2Kl)y^ - klx^) 

h = yil ~ (x^ + y^)) ~ xia - (x^ + i/)) 



(12a) 
(12b) 

(13a) 
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FIG. 5: Hopf-system with strong stochastic contribution along the j/-direction. Drift functions are (a) D[^\x,y) and (b) D2^\x,y). (c) 
Probability density function of {x, y) shows the most visited regions in phase space. Diffusion functions are (d) Jj^j' {x,y), (e) D['^ {x,y) — 
D^f {x, y) and (f) 1)':^^ [x, y). Results were obtained by analyzing the time series x = r cos 9 and y = r sin 6 according to Eq. Jilt with the 
same upward shifting for clarity as in Fig.|4] The probability density function in (c) explains the deviations observed for the diffusion functions 
(see text). 



+ ((^2 - 2Kf)x' - kly^) . (13b) 

Here we chose Ki ~ 0.05 and K2 — 0.5. In other words, 
the impact of stochastic fluctuations is large in the y-direction 
and small along the x-axis. Notice that the determinant of 
the diffusion matrix is not preserved, since the Jacobian of the 
transformation is not the identity matrix (see Append. |A). 

As already done in the previous example, we integrate sys- 
tem ( fTTT ) numerically and analyze the resulting data series of 
X and y values. As can be seen from Fig.|5] the drift functions 
hi and h2 are properly derived (see Figs. |5j and|5j)), as well 
as the terms of the diffusion matrices _D^^^ (Figs. |5}l-f)- The 
deviations observed for the small at particular regions of 
phase space (Fig.|5}l) are due to the lack of accurate statistics 
at those regions (see the PDF in Fig.|5};)- 

Interestingly, in both cases presented above, our analysis 
estimates better the diffusion matrix than the corresponding 
drift vector, contrary to what is known in other situations 
The reason for the poor quality of the drift estimate is probably 
due to the small time increment used for the estimate, which 
spoils particularly the result for the drift. See Ref. lfl9l for 
details. 

As can be seen in Fig. |6^ and|6j5, the eigendirections are 
properly derived as well as the two eigenvalues plotted in 
Fig.|6}; and|6}l as a function of r and 9, respectively. 

With these two examples we first have illustrated that 



our numerical approach succeeds in estimating the drift and 
the diffusion functions contributions in a two-dimensional 
stochastic system. A straightforward computation of the 
eigenvectors and eigenvalues of the diffusion matrices then 
enables us to ascertain a set of principal stochastic directions. 
We would like to note that the method described here bears 
some resemblance to the well-known principal component 
analysis i20il. In contrast to principal component analysis per- 
formed on the distribution of the measured data the method 
we propose, however, focuses on the stochastic fluctuations 
in time. The properties of these fluctuations in general are 
dependent on the position X resulting in vector fields of the 
principal axes of the stochastic contributions to the dynamics 
of the system of consideration. The two methods seem how- 
ever related and should provide complementary insights when 
applied to specific sets of empirical data. Such matters are be- 
yond the scope of this paper and will be addressed elsewhere. 

We also emphasize that the minimal eigenvalue of the dif- 
fusion matrix not necessarily points towards the direction the 
systems as a whole can be expected to evolve to. It rather re- 
flects the direction of minimal stochastic forcing only. For 
prediction of the system's evolution in general time propa- 
gators need to be taken into account which involve both the 
stochastic (diffusion) and the deterministic (drift) parts of the 
dynamics ||2|l- Nevertheless, being able to ascertain the direc- 
tions in phase space where fluctuations are weaker enables one 
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FIG. 6: Results of the analysis of diffusion matrices plotted in Fig. [5] obtained from a system witii an increased level of dynamical noise on 
the y component (second example). The individual panels depict (a) the vector fields of eigenvectors for the minimum and (b) the maximum 
eigenvalue, normalized to the respective eigenvalues, and the dependence of the eigenvalues both both on (c) r and (d) 6. 



to choose a transformation such that part of the new variables 
have small stochastic terms, reducing the number of variables 
which are affected by stochastic forces. Further, this is also of 
physical interest to see in which variables noise is acting, for 
example to which variable a thermal bath is physically con- 
nected, and it may be of technical interest to better detect a 
noise source, for example in an electric circuit ll2T1l . 

IV. DISCUSSION AND CONCLUSIONS 

In this paper we introduce the concept of eigendirections 
for the stochastic dynamics in systems of arbitrary dimension. 
The procedure builds on the modeling of complex systems 
by means of drift an diffusion functions, that specify a sys- 
tem of coupled Langevin equations. Estimates for these func- 
tions can be obtained directly from measurements on the sys- 



tems without prior knowledge on its dynamics following an 
approach described in Ref. |lll|2|l. 

In Langevin systems the stochastic forcing is composed 
of independent Gaussian i5-correlated stochastic fluctuations. 
The way how these fluctuations effect the system depends 
both on the diffusion matrix and on the coordinate system cho- 
sen. 

Whereas in former publications much attention has been 
paid to imperfections of the Langevin process, like measure- 
ment noise or correlated noise and finite size sampling ||22l - l25ll 
here we introduced a method which allows us to determine the 
eigendirections along which each stochastic force acts. Each 
direction of forcing is defined through the eigenvector and the 
corresponding eigenvalue accounting for its amplitude. The 
set of eigenvectors does not depend on the choice of the co- 
ordinate system and is therefore characteristic for the system. 
Further, for the particular case where a number k of eigenval- 
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ues are negligible in comparison with the others at each mesh 
point, the number of stochastic variables can be reduced. Even 
in cases where the number of stochastic variables cannot be re- 
duced, the eigenvector associated with the lowest eigenvalues 
at each point in phase space indicates the path with minimal 
stochastic forcing. In any case a transform to new coordinates 
in the directions of minimum stochastic forcing can be per- 
formed, increasing the relative amplitude of the deterministic 
components and the predictability of the corresponding vari- 
able. 

The method was successfully applied to a two-dimensional 
system exhibiting a Hopf-bifurcation. For the Hopf- 
bifurcation the diffusion matrix is better estimated through our 
analysis than the corresponding drift vector, contrary to what 
is known in many other situations 0]. 

Compared to previous methods for minimization of 
stochasticity introduced in Refs. JHIl^, our approach has the 
advantage of not requiring a parametrized Ansatz for quan- 
tifying the respective stochastic contributions to the system. 
Moreover, it should be applicable even in the case where the 
data sets are contaminated with measurement noise 
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Appendix A: The diffusion matrix in an arbitrary coordinate 
system 

Consider a transformation of variables X = {Xi}— >X = 
{Xi} with i = 1, . . . , A^, which is given by a two-times con- 
tinuously differentiable deterministic vector function F 

X = F(X,t) . (Al) 

Using Ito's formula lfl2l[l3ll . a truncated form of the Ito Taylor 
expansion, the Langevin equations for the new variables take 
the form 

dX, = ^dt + J2^dX,+ 
ot dXk 

^ N N g2p 

k=l 1=1 

Eq. (HI can be rewritten to 

N 

dX, = Df\l^)dt-vY.g^JOqTJ{t)Vdt . (A3) 



Inserting this expression into Eq. (IA2) . retaining all terms in the expansion up to order dt while neglecting terms of order 
{dtY'/'^ or higher, and taking advantage of the statistics of the stochastic forcing, Eq. (|5]l, one obtains for the third expression on 
the rh.s. 



JVAf TV „2p ^ N N N 

J„ 1 1 1 A 1 ^ ^ ™ 1 J 1 ( 1 A ™ 1 



Thus, the transformed Langevin equation has the form 



k—l l — l j,m—l 



'dXkdX, 



■5jradt (A4) 



^ A 1 I 1( I.-,- 1 ^ ^ 



fc=l 



k=l 1=1 j=l 



(A5) 



Restricting to stationary processes ( = 0) and using the The components of matrix G transform as 



notation 



dJQ 
dt 



N 



7)f^(X) + E.g,,(X)r,(t) (A6) 



the transformed drift function reads 



fe=i 



gijgkj 



1=1 3 = 1 



dXkdXi 



(A7) 



N 



9ij = E ^''^djT ^ G = JG 



k=l 



(A8) 



where J is the Jacobian of our transformation from coordi- 
nates Xi (vector basis e^) to Xi (vector basis e^). 

From elementary algebra it is known that a vector v can be 
written in both coordinate systems: 



^ = J2 "^'j'^j = E E 



k^k 
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(A9) 



A; 

E 



Thus, the matrix U incorporating the columns of eigenvec- 
tors Uk of matrix G, with coordinates in basis e^, can be writ- 
ten as U = J^U with U = [ui U2 ... un] fulfilling 
= U^. Thus, from Eq. ( IA8b one obtains 



U'GG'U = U'GG'U. 



(AlO) 



So, the eigenvectors of D^^' and D^^' are the same, apart 
from the basis in which they are considered, and in the orig- 
inal (Cartesian) coordinate system X, the transformation to 
principal axes is given by Eq. (lAlOl i i.e. 



U^D(2)u = diag[Ai A2 



A 



(All) 



At regular points the transformation in Eq. (lAU can be in- 
verted, 



X = F~\X) =: f(X) 



(A12) 



By definition f (X) is chosen such that the normalized eigen- 
vectors are given by 



Ufc 



1 df 



with Sk ~ 



df 



dXk 



sk dXk 

The Jacobian of f (X) can then be written as 

or in a more convenient form 
J = [ui U2 ... un] diag[si S2 ■ ■ 



(A13) 



(A14) 



sn] —'■ Us. 
(A15) 



The diagonal matrix s describes the metric of the new sys- 
tem. The Jacobian of the transformation in Eq. (lAlt is the 
inverse relation 



J = s-iU^. 
Introducing this relation into Eq. (IA8b yields 
G = JG = s^U^G , 
i.e. the transformed diffusion matrix reads 

D(2) = GG^ = s-^U^GG^Us-i 



(A16) 



(A17) 



(A18) 



Inserting Eqs. (lAl lb and (IA15b finally leads to the diagonal 
matrix 



D(2) = diag 



Ai A2 



L''l 



(A19) 




FIG. 7: Interchanging between coordinate systems. Arrows indi- 
cate that eigenvectors of the eigenvalues of the diffusion matrix (see 
Eq. l |A19b ) at each point of the mesh from which drift and diffusion 
coefficients are derived. Here the two-dimensional case is illustrated. 
Extension to arbitrary number of variables is straightforward. 



In practice, we have no access to the transformation in 
Eq. (lAlb . Instead we have a grid (numerical) representation 
in X coordinate system and search for the transformation in 
Eq. (lAlb is such that one of the new variables, say Xi, is 
defined by the field of eigenvectors Ui corresponding to the 
maximal eigenvalue Ai. Then, each orthogonal direction is 
consequently defined by one of the other (orthogonal) eigen- 
vectors Uj with j ~ 1, . . . , N. Figure |7]illustrates this for the 
two-dimensional case. 

To obtain the coordinates of each point P in this grid in the 
new coordinates X one considers the transformation (lAlb to- 
gether with the inverted transformation (IA12t and solves the 
system of PDEs in ( IA13b . The problem lies in the coupling 
introduced by the factors Sk and the additional complication 
of finding the correct boundary conditions which turns the so- 
lution of the PDEs into a complicated problem, even numeri- 
cally, although in 2D the can be neglected as we are only 
interested in finding the direction -§^- In most cases, how- 
ever, it should be possible to guess a suitable transform from 
visual inspection of the fields of eigenvectors. From one point 
to the next one in the mesh, the eigenvectors are sorted ac- 
cording to continuity arguments. Further, depending on the 
obtained fields of eigenvectors, a scaled polar form or hyper- 
bolic coordinates may be considered, where the parameters of 
said transform can then be fitted to the vector fields. Such par- 
ticular cases depend on the specific data set at hand and will 
be addressed elsewhere. 
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